Sedimentation and Flow Through Porous Media: Simulating 
Dynamically Coupled Discrete and Continuum Phases 



Stefan Schwarzer* 

Laboratoire de Physique Mecanique des Milieux Heterogenes, Ecole Superieure de Physique et 

Chimie Industrielles, 
75231 Paris, Cedex 05, France 
Hochstleistungsrechenzentrum, Forschungszentrum Jiilich, 
524-25 Jiilich, Germany 
(February 1, 2008) 



Abstract 



We describe a method to address efficiently problems of two-phase flow 
in the regime of low particle Reynolds number and negligible Brownian mo- 
tion. One of the phases is an incompressible continuous fluid and the other 
a discrete particulate phase which we simulate by following the motion of 
single particles. Interactions between the phases are taken into account us- 
ing locally defined drag forces. We apply our method to the problem of flow 
through random media at high porosity where we find good agreement to 
theoretical expectations for the functional dependence of the pressure drop 
on the solid volume fraction. We undertake further validations on systems 
undergoing gravity induced sedimentation. 
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I. INTRODUCTION 



A classical problem of chemical engineering is the understanding of particulate two- 
phase flows, in which a continuous fluid constitutes one component and a discrete particle 
phase the other. The practical importance of understanding such a particle laden flow is 
evidenced by its central role in geophysical phenomena like sand storms, dune formation, 
or sediment transport, by the importance of biological questions like the understanding of 
the flow properties of blood, or cell component separation techniques as ultracentrifugation, 
by industrial applications as diverse as pneumatic transport in tubes, catalytic cracking, 
biological and chemical reactors, solid fuel rocket motors, fluidized beds, sedimentation, 
filtration and many more |],|2[]. 

If the discreteness of the particulate phase is fully taken into account then the mathe- 
matical formulation of the flow problem involves (i) a field equation for the continuous liquid 
phase — in most cases the Navier-Stokes equation, subject to a set of boundary conditions 
both on the container walls and the surfaces of the suspended particles — and (ii) a set 
of differential equations for the time evolution of the degrees of freedom of the individual 
constituents of the particulate phase. Due to the extremly complicated boundary conditions 
and the nonlinearity of the underlying equations, an analytical solution to this problem is 
impossible, except for exceptionally simple cases. A numerical treatment, however, even 
with simplifying assumptions, still poses tremendous practical problems. 

Several simulation techniques have been developed which we briefly review here, (a) 
Finite volume techniques that implement no-slip boundary conditions on the surface of each 
particle have been employed for very few particles. The treated Reynolds numbers |§ are in 
principle only limited by the grid resolution and the computation time available. Similarly, 
(b) lattice Boltzmann techniques for the liquid equation have been used to simulate up to 
1024 suspended particles Q. So far, these have come closest to provide realistic simulations 
of two-phase flows with many particles. Most other approaches involve certain assumptions 
that are only true in limited parameter ranges: One important class of algorithms uses the 
Stokesian dynamics technique [|||| valid for low Reynolds number flow to (c) obtain the flow 
field consistent with no-slip boundary conditions J7| or to (d) approximate the particles as 
being point like 

(e) More popular in the engineering sciences, but less rigorous are continuum approxima- 
tions that involve two sets of continuum equations, one for the liquid phase and one for the 
particulate phase. The boundary conditions on the particle surface and their influence on 
the flow is then represented by local drag forces depending on the solid volume fraction and 
local velocity differences [p~0| , |TTj| . In these approaches there are remaining open questions 
in the determination of the proper constitutive equation for the solid and the momentum 
exchange between the phases. 

(f) To circumvent some of these later problems of the pure continuum approaches, algo- 
rithms have come into fashion that combine a discrete element |12[ or molecular dynamics 
description [T3|] of the particles with a continuum equation for the liquid as in (a), but use 



a drag term to couple the two phases as in the pure continuum formulations (e) [[14| ,|T5[ 

In this paper we will describe an algorithm of type (f). Such a method allows imme- 
diate access to basically all physically relevant quantities in the system, including particle 
coordinates and both particle and liquid velocities, at computational costs comparable to a 
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"standard" real space Wavier-Stokes integration. The main drawback is that a neglect of 
the proper boundary conditions in the treatment of the liquid will result in an inaccurate 
rendering of the short scale flow properties. Since, however, our main focus will be the 
ability of the algorithm to describe collective phenomena, i.e., the effects that arise when 
the number of particles is large, we do not have the ambition to describe accurately the local 
flow fields on the scale of the particle size. At the moment the detailed effects of neglecting 
certain aspects of the local flow are not clear to us. Nevertheless we want to see which and 
how some collective phenomena emerge directly from simple modeling assumptions — as 
opposed to using semi-empirical expressions as, e.g., done in Refs. |T3,|T3[. In particular, we 



rely on the fact that the long-range hydrodynamic interactions, which we presume to be the 
most important for collective phenomena, are correctly represented by the velocity and pres- 
sure fields according to the Navier-Stokes integration under consideration of an additional 
field representing the particle density. 

The purpose of this paper is to give a detailed description of the simulation algorithm 
that we use to model particulate two-phase flow and its validation in the cases of (i) flow 
through a "porous medium" consisting of the "particles" in the simulation which are kept 
at fixed positions and (ii) of particles that undergo sedimentation in a suspension under the 
influence of gravity. In Section [TT] we will first describe a two-dimensional implementation 
of the algorithm. Then, in Sec. [TTI] , we discuss first the case (i) of flow through a porous 
medium (Sec. [Ill A|) and second the case (ii) of sedimenting particles (Sec. [IHB|). 



II. THE ALGORITHM 

We aim at a simulation of macroscopic, hard sphere particles that interact on contact 
when the interstitial fluid plays no important role. In addition, they interact over long-range 
hydrodynamic forces, mediated by an incompressible fluid medium between the particles. 
For the moment, we consider a two-dimensional implementation of our project. Since the 
analytic aspects of the solution of the Navier-Stokes equation in two dimensions are very 
different from the three dimensional solution, we expect at best a qualitative agreement 
to the features of real three dimensional flows. However, we feel justified to study the 
two-dimensional problem if it is possible to find qualitative knowledge on that way. 

We organize the description of our algorithm into three main parts. In the first section 
we describe details of the employed molecular dynamics technique for the particulate phase. 
The second compiles the fundamental equations for the liquid phase and the details of their 
solution. In the third section we elaborate on the particle-liquid interaction. 



A. Molecular Dynamics of the Particle Phase 

1. Single particle properties and forces 

Each particle i in the simulation is characterized by its radius r^. We take the particle size 
distribution to be slightly polydisperse with the radii r» drawn from a Gaussian distribution 
h p (r/r) with mean f, which is cut off at its standard deviation pf, 
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h p (r/f) oc < L 2 J ' '''' (1) 



Here, the distribution is written in terms of the dimensionless radius r/f and the polydis- 
persity parameter p. 

We disregard effects of rotational motion of the particles and only consider their trans- 
lation. Thus we describe the particle motion only by the coordinates of their center of 
mass Xi = (xi,yi). Although the particle centers are constrained to lie in the xy plane, we 
will assume — for reasons that we address in Sec. [11 C| — that the particles can other- 



wise be regarded as three dimensional. Therefore we associate with each particle a mass 
rrii = (4/3)iirfpp, where p p is the constant particle density. 
We consider the following forces to be acting on particle i, 

Fi = \^{ Pp - Pl )g + F? + F? + + F*). (2) 

Here, Fi is the sum of all forces on i, the F™ are forces due to the presence of the boundary, 
the term proportional to g accounts for particle weight and buoyancy. The drag force Ff 
will be discussed in Sec. II C. The sum runs over the other particles in the system, but 



is effectively restricted to the neighborhood of particle i, since we will only introduce short 
range interparticle forces in radial (F£j) and and tangential (F^) direction in the next two 
sections. 

Container walls, if present, manifest themselves by elastic forces and frictional forces in 
analogy to those acting between two particles. We will use for them the equations that 
result from those for two collision partners — to be introduced below — in the limit of the 
other particle having infinite mass and radius. 

In direction of g, we account for the weight of the particles. Acting opposite to gravity, 
we add the buoyancy reaction of the liquid, which equals the weight of the displaced fluid 
of constant density p\. Our coordinate system is chosen such that g points in —y direction. 

The forces from Eq. (Q) enter the equations of motion for the particulate phase, rriiXi = 
Fi. We solve this coupled system employing a fourth order predictor-corrector algorithm as 
described, e.g., in [ |l"3|j . 



2. Pairwise interparticle forces in radial direction 

At very low Reynolds numbers the lubrication forces between particles play an important 
role. Since they have a divergent character when the particles come close, particles in 
principle do not touch in this regime. However, if the fluid in the system is a gas or the 
particle based Reynolds numbers are not any longer small or if the particles are assumed to 
have some surface roughness, then they will touch and we have to model the contact forces 
between them. 

If no fluid were present then the particles are force free unless they touch, whereupon 
strong repulsive and dissipative forces act between them, resulting from the viscoelastic 
properties of the particles. To model these forces, we here use contact models frequently 
employed in granular matter research P2]JT^JT7[ ]. The first of the forces acting on member % 
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of a particle pair i, j in contact is an elastic restoring force F-j. This force is proportional 
to the virtual overlap £y = (r* + rj) — \x{ — Xj\ of the particles. If the overlap is positive, 
then 

F$ = -kn^jUij, for > 0. (3) 

Here, = (xj — Xi)/\xj — Xi\ denotes the unit vector pointing from the center of particle 
i to the center of particle j and k n is the stiffness of the contact, which we assume to be 
constant.[] As mentioned above the force vanishes if the particles do not overlap, i.e., £y < 0. 

To take the dissipative character of the contact into account, we add a velocity dependent 
friction term to the elastic restoring force. This damping force shall also act in the direction 
of the line connecting the particle centers (the normal direction) and be proportional to the 
normal relative particle velocity. On inclusion of this term, we obtain for the contact force 
Fff 1 in normal direction, 

F%j ( ^n£ij Tn^red (Xj Xi) • riij)riij. (4) 

In this relation, m re d is the reduced mass miW,j/{mi + rrij) of the pair in contact and j n 
determines the strength of the dissipation. We have suppressed the indices on m re d and j n , 
but we consider them to vary among different particle pairs. 

Equation (f|) is the equation for a damped harmonic oscillator while the particles are in 
contact. For a given initial normal relative particle velocity it can be solved analytically 
for the velocity after the contact vj. Since energy is dissipated, the ratio of these velocities 
is less than one; its value is termed restitution coefficient e = I'UjfI/I"^!. One obtains for a 
specific particle pair (pair indices suppressed), 



exp 



± A. _ , 

V In m red 



(5) 



However, in the presence of a liquid at low Reynolds numbers, the dominant force at 
small distances between pairs of particles is the lubrication force F 1 arising from the pres- 
sure necessary to replace the liquid within the gap between the two approaching spheres. 
The lubrication force damps the relative motion of the particles very strongly and diverges 
when the particles touch. Since in our approximate treatment of the liquid flow in partic- 
ular the lubrication force is not well represented we add it as an additional component of 
the particle-particle interactions, only active at short distances between the particles. The 



normal component of the lubrication force is |18[j 



The determination of k for realistic contacts in 3D is not a trivial matter. The value of k 
depends on the exact physical processes governing the contact and thus in general both on material 
constants as the Young modulus E p or the Poisson number cr p and the geometry of the particles. 

~* 1 3/2 -* 

Many researchers use the nonlinear Hertzian contact theory, in which Ff- = —k n ^ Uy. For equal 
sized spheres, the Hertzian theory gives k n = V2rE p /3(l — ai). We refer the reader to the tables 
| and [0| for values used in the simulations. 
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ph n = -Qnri-^-Kxi - Xj) ■ n i3 \n iy (6) 

In the above equation, 77 denotes the shear viscosity of the liquid and r re d the reduced radius 
r i r j/{ r i + r j) of the pair. Again we have suppressed the indices i and j for simplicity. The 
expression is only valid for small positive separations between the surfaces of the two 
involved particles. These are reflected in negative overlaps ^ of small modulus. 

To take into account the limited range of validity of Eq. (||), we cut off where 
> r Ted . To avoid a discontinuity in the force law, we subtract in Eq. (|6]) a constant 
equal to the value of the force just at this cutoff distance. 

Furthermore, we remove the divergence at £y = by adding to the value of (— £y) in 
the denominator a small positive number 5r red . We have used a value of 5 = 0.1. Other 
than just being a numerical contrivance our physical motivation is found in the unavoidable 
surface roughness of particles in reality, which may cause the particle to come into contact 
despite the lubrication force. 

Similarly, contacts due to numerical inaccuracies are, almost unavoidable in a dense 
system with many particles. To cover these spurious cases as graceful as possible, we have 
matched the 7„ of Eq. such that the force law is continuous when £ = 0, for particles of 
radius f. We obtain, 

7n = Qtt V ^- (- - j^—) ■ (7) 
m red \0 1 + 0/ 

The verbal statements in the preceding three paragraphs condense into the following 
equations for the total interparticle forces in the normal direction, comprising both contact 
and lubrication forces, 

fin = fic,n fil,n 
tj l l 3 

{ , if > r red , 

= | [-^VrlM - Si) ■ n l3 ][ _^l SrrcA - {1+ s )ri J n ij: if < (-&) < r red , (8) 

[ [ — kn^,ij ~ ln m red(%j ~ x i) ' n ij)] n ijy if £ij > 0. 

This force is continuous over the whole range of £y. 



3. Pairwise interparticle forces in tangential direction 

To model the frictional forces acting perpendicular to the line connecting the two par- 
ticle centers — the tangential forces — we resort, as in the case of the normal forces, to 
notions of particle contact modeling. Here, the Coulomb law of sliding friction asserts that 
the magnitude of the tangential friction force F\- is — on contact — proportional to the 

— # 

magnitude of the acting normal force F£j, [from Eq. (§)] with a constant of proportionality 
\x usually between 0.05 and 0.5, 

l4i=^% ( 9 ) 

This force is always directed opposite to the relative motion. Numerical problems may occur 
in near central impact, when the tangential component of the relative velocity is small, but 
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— * 

F n is large. Then the likewise large tangential component of the force resulting from (£J) 
may cause an unphysical oscillatory behavior of the tangential velocity during contact. We 
therefore replace Eq. ([|) for small relative tangential velocities by a velocity proportional 
friction term. Thus, finally, the tangential friction force on particle i becomes 

->t 

4 = -min(/i|it|, tH'DlSt' ( 10 ) 

\ v ij\ 

where ■ has been introduced as an abbreviation for the relative tangential velocity (xi — 
Xj) — [(Si — Xj) ■ nij)fiij. We do not include shear contributions of the lubrication force into 
the interparticle forces. For simplicity, j t is taken to be constant. 



B. Fluid Model 



We describe the state of the fluid phase by three continuum fields, namely (i) the velocity 
field u(x) of the fluid, (ii) its pressure p(x) and (iii) a field e(x) equal to the local volume 
fraction of liquid. These variables have only physical meaning as averages over volume on 
a scale larger than that of the individual particles. Their choice is motivated by continuum 
approaches to multiphase flow [ |T0|| . 

The position and geometry of the particles determines the field e(x) which — for specific 
discrete tiling of the simulation plane — is a more or less smooth function varying between 1 
(no particles) and (full occupation by particles), defined for all tiles and all times. Similarly, 
the time evolution of the velocity field u(x) is determined by the pressure distribution, 
viscous contributions and a force distribution f(x) which comprises both volume forces on 
the liquid and momentum exchange contributions with the particulate phase (Sec. [11 (J| ). 



We follow Ref. |10| and write for the time evolution of the liquid velocity u (we will drop 
the argument x of the fields from now on) , 



- + (U V) U 



-eVp + er] V 2 m + ef. 



(11) 



Although e drops out from this equation, it enters into the the momentum exchange contri- 
bution to /, in the sense that the momentum transfer to the particulate phase due to the 
fluid phase — due to drag between particles and liquid — must be "fed back" to a liquid 
volume smaller than in the case without particles. 

Since we have in mind applications to systems with typical velocities much smaller than 
the velocity of sound, we can assume that the fluid phase is incompressible. The equation 
of liquid mass continuity then reads 



depi 
dt 



+ V 



eu 



0. 



(12) 



Equation (|T2| ) presents a constraint on the velocity field that must be fulfilled at all times 
and may be employed to obtain the pressure field via an iterative procedure which we model 
after the artificial compressibility method of Chorin [|T9|-f2T|j . Here, we sketch the basic ideas 



of its 2D implementation briefly and refer the reader for more details to the literature. 
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We discretize the differential equations ([11]) and (|T^) in the following way. We place the 
velocity components well as the pressure p on three quadrilateral meshes with 

lattice spacing Ax. With respect to the pressure grid, the grids for the x and y velocity 
components are shifted by Ax/2 in x and y direction respectively. This construction is 
commonly refered to as the MAC mesh and has several computational advantages. For 



instance, it is a simple means to avoid numerical instabilities due to mesh decoupling |2l 
The choice of location for the computational quantities is conceptually related to location 
of variables in finite volume techniques for flux conservative differential equations, in which 
the fluxes are located on the corresponding faces of a control volume whereas the conserved 
quantities themselves reside in the center of the volume . 

We obtain the pressure and the velocity components by an iterative procedure. Let the 
index n refer to values at time t = t n and n + 1 to those at t = t n+ i = t n + At after a timestep 
of duration At. The index k shall denote an iteration index. We define p n +i,o = Pn, i-e., 
we start an iteration for the new pressure at time t n+ i with the old values at t n . We obtain 
a tentative velocity field at t = t n+ i from an evaluation of the discretized Navier-Stokes 
equation (|TT|), 

Pi — = -P\{UrN)u n ~ Vp„,+l,fc + f n + r]V 2 u n . (13) 

where the symbol V now denotes second order precise difference operators on the lattice. 
As mentioned above, the solid volume fraction enters implicitly through f n . However, since 



(|13|) is a discretized Navier-Stokes equation, its stability criteria on the MAC grid have been 
studied and are reported, e.g., in These criteria state that the values Ax and At are 

subject to the two constraints 

At < -7i P~i ivr> ( 14 ) 

~ Pi(|< ax | + |w™ ax |) 2 



and 



At < (15) 
4?7 



In general, the velocity field u n+ i t k+i resulting from ([13]) considered together with e n does 
not satisfy the continuity equation (|12"D. Rather, one has to conceive the continuity equation 
as a constraint that determines the pressure field within the fluid such that the resulting 
velocity field satisfies the continuity equation at all times. 

To this end, one derives an iterative procedure to determine an appropriate pressure 
field. This procedure is based on the idea that the local violation of the continuity equation, 
i.e., the value of the left hand side of Eq. (|T^), can be used to correct the pressure field. 
The correction is taken in a direction such that the modulus of the violation is reduced in 
the next iteration step, after a new tentative velocity field has been determined. We write 

Pn+l,k+l = Pn+l,k - Vl[-7^- + V • (e n U n+hh+ x)}. (16) 

Here, a large value of the parameter A is crucial for rapid convergence. The value of A is, 
however, by stability requirements constrained to 



S 



(17) 



In the simulations we have chosen A to equal its upper stability limit. 

It should be noted that the values e n in Eq. fll6|) are not all located on the same subgrid: 
the time derivative is taken at the pressure points and the e n multiplying -u n+l fc+1 is realized 
by different fields, each living on the same subgrid as the associated velocity component. 

Once we have determined a new pressure field, we need to recalculate new velocity 
values consistent with the p n+ i^+i- These velocities may be obtained using the Navier- 
Stokes equation (|I3"D , or equivalently, avoiding the costly reevaluation of Eq. (|13|), using the 
relation 

— * — * 

Pi ^ = V(p n+ i,fc+i - Vn+l,k)- (18) 

Iterating Eqs. (|IED and ( |T8"D yields a new pressure field and after convergence a velocity 
field for time t n+ i such that the equation of continuity is satisfied. 

For our purposes, the described algorithm has three advantages. It (i) generalizes 
straightforwardly to 3D and (ii) unlike spectral or streamfunction methods it gives im- 
mediate access to the quantities p and u in real space. Fast access to the latter is crucial for 
the calculation of the particle-liquid interaction. Moreover (iii), only the chosen coarseness 
of the spatial and temporal discretization limits the range of Reynolds numbers address- 
able in the simulation. However, since the presence of the particles is communicated to 
the liquid among others through the field e, it does not make much physical sense to use a 
computational grid on a scale smaller than the particle size. 

For Reynolds numbers requiring such smaller grids, one could think maybe of an ap- 
proach to decompose the liquid equation into an equation for the average flow and additional 



equations describing the fluctuations around it, in the spirit of turbulence modeling [24 . 



C. Interaction of Particles and Fluid 

The main problem in simulations of multiphase flow that try to bypass the specification 
of boundary conditions on the phase boundaries is to specify an expression for the mo- 
mentum exchange between fluid and particulate phase. Even for a single sphere this task is 
daunting [2l| . For extended fixed random collections of spheres the phenomenological Ergun 
formula [PIP?! gives the pressure gradient as a function of solid fraction and liquid velocity. 



Tsuji and coworkers fl4l , |T5|| have used a "localized" form of the Ergun equation to estimate 



the force of the liquid on a sphere and obtain for pneumatic transport in pipes and fluidized 
beds qualitative agreement of flow patterns and quantitative agreement of some quantities. 

Having the complicated situation in mind, and being aware of the significant approxima- 
tion involved, we wish to explore the consequences of a very simple model for the momentum 
exchange. We use a local version of the "global" Stokes formula — which is just one term in 
the expression given by Maxey for the drag on an isolated sphere under low Reynolds 
number conditions — to obtain the drag force acting on a sphere of radius r^, 

F? = -during - u(xi)}. (19) 
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In order to evaluate the liquid velocity field at the location of the particle we interpolate 
linearly the velocity values from the closest four grid points. That is to say, if the pair 
(x, y) denotes the coordinates of one of these four "adjacent" grid points Xj then w(x,y) = 
(1 — \xi — x\/Ax)(l — \yi — y\/Ax) is the weight associated with (x, y); the index i refers to 
the particle location. We obtain thus for some velocity component u the interpolated value 
u(xi) as u(xi) = ^w{x,y)u{x,y). The sum extends over the four corners of the MAC grid 
plaquette associated with the component u into which Xi falls. 



As equation (|T^) resembles the 3D expression for the drag on an isolated sphere, we should 
give here a motivation for this modeling assumption. A "pure" 2D simulation should, strictly 
speaking, be one of rigid, parallel cylinders. However, at a fixed small Reynolds number, the 
drag per unit length of a cylinder does not depend on its radius. This behavior is very 
different from the observations in 3D experiments, where there is a strong dependence of 
the drag on the size of the spheres. Thus, in particular, effects of the particle polydispersity 
cannot be expected to be represented properly by a model based on drag-forces of cylinders. 

We therefore use Eq. (|Tj|) for the drag and refer the force ([19]) to a reference length z 
equal to the average particle diameter. The picture we have in mind behind this procedure 
is that we imagine the particle configuration repeating itself all z length units, and the liquid 
being infinitely extended in z direction, however both particles and liquid being constrained 
to motion in the xy plane, when we calculate the drag per unit volume that enters into the 
Navier-Stokes equation. The finite drag force ( |I"9D acts on the liquid at the location of the 
center of the particle and we must thus — in order to convert it to the force density required 
in the continuum equation — represent it by a 5- function at Xi, 

Ff^—Six-Xt), (20) 
ze(x) 

in the Navier-Stokes equation. The e(x) ensures that the force density is refereed to the 
liquid fraction alone which is necessary to conserve momentum. We form the sum over all 
particles and add the uniform contribution of gravitation to obtain the full volume force 
density term, 

fXx) = g + Y,F l d ^—Mx-x l ). (21) 



On our computational lattice we implement the expression ( |2"T| ) by distributing 
Ff I z(Ax) 2 e(x) to the four grid points closest to Xi. To this end, we employ the same 
weights as introduced for the interpolation of the velocity components above. For example, 
the contribution to point (x,y) is (1 — \xi — x\/Ax)(l — \y.- t — y\/Ax)Ff/z(Ax) 2 e(x). 

This concludes the description of the implementation of our algorithm. 



III. APPLICATION TO EXAMPLE PROBLEMS 

We will now describe the results of the application of the described algorithm on selected 
physical problems in order to validate our approach, and assess its limitations or respectively 
the consequences of using the simple drag law (^9|). 
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A. Flow through porous media 



As our first system, we consider the motion of a fluid through a fixed random assembly 
of particles. Our mental picture is that the particle assembly is a model for a random porous 
medium with very high porosity. In the simulation, we then have to keep the particles pinned 
to their initial positions. However, we calculate all forces acting between particles and liquid, 
but we only update the liquid's degrees of freedom in the computational timestep. Excluded 
volume effects due the presence of the particle phase and local frictional drag still influence 
the fluid motion. Gravity is set to zero, or equivalently the flow is considered to lie in the 
horizontal plane. 

In x direction the system has the width L x and boundary conditions are periodic. We 
impose fixed superficial flow velocities u in y-direction at the inlet and outlet of the system 
at y = and y = L y , where L y denotes the height of the system. A typical arrangement 
of particles in a small system of L x /f ~ 33 and L y /L x = 2 at Re p = 3.75 x and the 
resulting stationary flow pattern is displayed in Fig. [I]. The circles indicate the particles, 
the arrows direction and magnitude of the liquid flux obtained by multiplying the local 
liquid volume fraction with the flow velocity. We see how the particle volume fraction 
influences the flow pattern such that the current concentrates in regions with few particles 
present. Note that the fluid velocity is defined everywhere in space, even at points covered 
by particles and that the flow velocities should therefore be considered as average values 
over the specific grid cell. Note also that the flux vectors are not displayed at their location 
used for computational purposes but are extrapolated to and displayed together with the 
color coded pressure at the location of the pressure points. 

In this system, we measure the overall pressure drop per length Ap/L y as a function of 
the overall liquid volume fraction e of the medium and the superficial liquid velocity u. In 
the viscous regime the pressure drop is proportional to the fluid velocity u. We evidence a 
constant ratio of pressure drop to fluid velocity in Fig. ||] for several orders of magnitude of 
the particle Reynolds number Re p = pifu/rj at fixed "porosity" e. 

We report our results in terms of the friction factor 

f p = -fAP/ Pl u 2 L y , (22) 

which is, due to the linearity of the pressure drop in the viscous regime, inversely proportional 
to the particle Reynolds number Re p = p\fu/rj. Accordingly we have plotted in Fig. |3] the 
product f p Re p which — in the viscous regime — is independent of the Reynolds number. 
The product's value is plotted against the liquid volume fraction or porosity e based on the 
3D volume of the particles divided by the box volume L x x L y x z, where z denotes the 
effective box depth introduced in Sec. [II Q . The three different curves in the plot indicated 



by symbols correspond to two different values of z/f = 2 and 8/3 and one run (lowest lying 
curve) where only drag between particles and fluid was considered in the simulation and the 
local liquid fraction e(x) in Eqs. ( |TTD and ( |T2"D was kept equal to 1. The curves for different 
z/f collapse into a single universal curve. We consider this data collapse as an a posteriori 
justification of the picture of the simulated system which we have used to find a form for 
the momentum feedback to the liquid (Sec. [11 C| ). 

Concerning the functional form of the result in 3D, the literature lists several phenomeno- 
logical expressions for the porosity dependence of the friction factor |29||. Popular expressions 
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for f p Re p include (i) the half-empirical Carman-Kozeny relation f p Re p ~ (1 — e) 2 /e 3 and (ii) 
the phenomenological Rumpf-Gupte form f p Re p ~ e~ 55 . Both relations are only valid in an 
intermediate range of porosity and do not apply to the limit of very high porosity or very 
low solid volume fraction, where one expects f p Re p ~ 1 — e to leading order for the following 
reasons. For a strongly "dilute" random medium, let us consider the solid as independent 
spherical "defects." Each defect exerts the Stokes force 67rrr]u on the fluid, where u is the 
interstitial fluid velocity. Thus we approximate, on the one hand, the internal rate of energy 
dissipation to be 67rrr)u 2 N, where we denote the total number of 'defects' in the system by 
N. On the other hand, externally, the pressure drop across the system and the superficial 
fluid velocity u = ue give the rate at which work is done on the system to be APL x zu. 
Equating the two rates and using (1 — e) = (4/3)7rr 3 iV/ \L x L y z) yields 

f P Re P = (23) 

This equation correctly predicts the pressure drop to vanish in the limit of large porosity, 
but it is clear from the argumentation given above that its validity is restricted to the regime 
of large liquid volume fraction. 

The solid line in Fig. |3] shows the prediction of formula (p3|). The simulation data 
compares very favorable with this prediction for the whole range of simulated effective 3D 
porosities e > 0.7. Please note that fl23"|) has no adjustable parameters. The agreement is 
surprising considering the rather crude assumptions that we have made. 

In contrast, the Carman-Kozeny formula predicts the pressure drop to vanish quadrat- 
ically as e — > 1, and the Rumpf-Gupte formula predicts even a constant pressure drop for 
e — > 1. However, both formulas are applicable only in the intermediate or low porosity 
regime so that the lack of agreement with (|23|) is not surprising. 

B. Sedimentation 

We now discuss the application of our algorithm to the more general case when both 
particles and liquid are mobile and important dynamical consequences arise from the cou- 
pling of the two phases. A prototype example for this case is batch settling sedimentation. 
Imagine a homogeneous mixture of particles and liquid being placed into a quadrilateral 
- or for our purposes rectangular — container and being initially at rest. If there exists 
a density difference between particles and liquid, then gravitational driving forces will set 
the particles in motion. Complicated fluid-mediated hydrodynamic interactions between the 
particles lead to convoluted trajectories and give rise to collective phenomena as for exam- 
ple anisotropic self diffusion of particles. The particles slowly settle to the bottom of the 
container with an average speed (V($)) that decreases as the volume fraction $ of particles 
in the container increases. We have chosen here $ = 1 — e to follow the convention in most 
of the literature on sedimentation. We display a typical situation during the batch settling 
process in Fig. [|. 

The conditions of the simulation are periodic or respectively no-slip boundary conditions 
in the horizontal x direction and no-slip boundary conditions in vertical ^/-direction. Gravity 
is taken to act in — y direction. One denotes as the hindered settling function /hs(3>) the 
ratio of the sedimentation velocity (V(<&)) to the Stokes velocity Vs, 
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(v)/v s = As($), 



(24) 



where 

V s =(2/9)(p p -p l )gf 2 /r ] , (25) 

which is the settling velocity of an isolated sphere in an infinitely extended fluid. 

We investigate /hs(3>) for an almost monodisperse system. The particle Reynolds number 
Re p = pifVs/t] is smaller than 1, typically ~ 3 x and (ii) the "simulation box" Reynolds 
number Re^ = p\L x Vs/r] is 100 times larger, i. e. Rey, ~ 3. Furthermore, we choose the 
particle size large enough such that the effects of Brownian motion may be neglected, cor- 
responding to the regime of high (iii) Peclet number. We have performed simulations and 
determined the settling velocity as a function of the solid fraction of the particle suspen- 
sion. The following four different sets of conditions for "thought" experiments have been 
performed to assess the role and importance of lubrication and backflow for the simulations. 
In particular, simulation series 

(i) includes lubrication effects [Eq. (^)], effects of particle void fraction [Eqs. ( P] , |I2] )], 
and periodic boundary conditions in x direction, perpendicular to gravity; 

(ii) differs from series (i) only in the respect that we have not considered the lubrication 
term (|6]). Without these lubrication forces, we choose the value of 7 n in Eq. (Q) such 
that for pair collisions a restitution coefficient of 0.9 results [cf. Eq. @]; 

(iii) differs from (ii) in the respect that we have additionally set the liquid fraction e(x) 
to 1, as if the particles consisted of liquid and not of a separate solid phase. The 
interaction of particle and liquid phase results only from the pointlike frictional drag 
between the two. 

(iv) differs from (i) only in the respect that we have used no-slip boundary conditions at 
x = and x = L x . 

Figure ^| displays the hindered settling function of case (i). The "Stokes" velocity used 
to normalize the data has been obtained from simulations of single spheres. In addition, we 
have plotted a theoretical expression (j3~0D, fhs(&) = (1 — $) 3 , which we will discuss later in 
the text (Sec. |1V Jj| ) . It should be noted that experiments in the viscous regime — both Re^ 



and Re p much smaller than 1 — often report a correlation with a (1 — $) n , n m 5 dependence 
3U| , |3"1"| | , the Richardson- Zaki formula [32|. It may be that in our case the larger value of 9ftb 



prevents us from seeing n « 5. Of course, also the nature of hydrodynamic interactions in 
2D is quite different from the 3D case, such that probably only full 3D simulations see n = 5 



The graphs in Figure |6| address the question in which way the resulting sedimentation 
velocities in series (ii) to (iv) differ from that in (i). To this end, we have calculated the 
ratios of the sedimentation velocities. First, the symbol (O) denotes the ratio (V^) / (V^) . 
Since lubrication alters the interactions between particles more significantly when the par- 
ticle fraction is high, we see that the fraction deviates more and more from 1 as the solid 
fraction increases. The ratio becomes smaller than 1 which means that the system "without" 
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lubrication sediments slower. We understand this behavior taking into account the strongly 
damping character of the lubrication force. Such Damping forces between particles have the 
effect to reduce the relative velocity between them and thus to favor the creation of loosely 
connected particle agglomerates. These effects are well known for "dry" granular systems 
with dissipative contact interactions between particles collisions Q,^]. These particle clus- 
ters then trap some liquid within and fall as almost coherent units. Such a unit displays a 
larger Stokes velocity than its constituent particles since Vs is proportional to the squared 
radius of the falling object. 

Second, we assess the effect of backflow in the simulation by setting e(x) = 1 in the 
Navier-Stokes and the continuity equation ( |H||T2| ). The particles are then effectively point- 
like as far as the fluid is concerned. Denoted by □, the ratio (V^ 11 ') / (V^) starts at a value 
below one — which results from particularities in the momentum exchange modeling (see 
Sec. [IV B| ) — and increases with volume fraction. We can estimate the effect of backflow 
by assuming that the average relative velocity of particles and liquid phase are the same 
in situations (ii) and (iii): A short calculation assuming that the average relative velocity 
between particles and liquid is the same in both cases predicts (V&*"))/(yM) = 1/(1-$) w 
1 + $. The observed increase is less steep as expected from this relation indicating a more 
complicated effect as consequence of the introduction of the volume fraction field, which we 
do not understand at this point. 

Third, we have performed a simulation with no-slip boundary conditions on the container 
walls instead of the periodic boundary conditions used otherwise. We denote the data points 
for (p'') / (yW) by the symbol O. The ratio is, apart from a point at very low volume 
fraction, smaller than one. The liquid velocity is constrained to be zero at the container 
walls and hence we believe that altogether the motion of the liquid is less vehement than in 
the periodic case. Since particle and liquid motion are very strongly coupled — the distance 
that a particle has to travel in order to reach a terminal velocity is much less than a particle 
diameter — the on average smaller particle velocities at the container walls suffice to slow 
down the settling. 



1. Relation of sedimentation and flow through a random medium 

At this point it is interesting to note that the functional forms proposed both for the 
volume fraction dependence of the mean settling velocity in the sedimentation problem and 
for the pressure drop in the problem of flow through random media (the Rumpf-Gupte 
relation) are probably not independent from each other. 

More precisely, it is possible to relate the formulas for the pressure drop in flow through 
a random medium and the settling velocity in the sedimentation problem. During the 
sedimentation process, one observes two shock fronts in the suspension. One front separates 
the densely packed particle assembly at the bottom of the vessel from the bulk of the 
suspension, the other the bulk from the leftover clear fluid region at the top. If the particles 
were frozen in their instantaneous positions by compensating the forces acting on them, then 
the pressure drop between these two shock fronts is just the hydrostatic pressure drop of 
the pure liquid. In an infinitely extended system, time averaged, each particle experiences 
a total force equal to the difference of its weight and buoyancy. If the particles move freely 
then an additional pressure difference between the two shock fronts arises because now the 
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liquid must support the particles. Therefore, in steady state, the arising pressure difference 
must just equal the difference of particle weight and buoyancy, i.e., 

AP 

— = $(p s - Pi)g. (26) 

Here L denotes the height difference within the two mentioned shock fronts over which the 
pressure drop is measured. 

On the other hand, if the particles were immobile, we could apply a drag formula of the 
form introduced in Sec. |III A| for flow through random media, say f p Re p = / rm (e), to find 
the pressure drop AP/L. Let in this case u denote the average relative velocity of liquid and 
particles. Then we find that 

— = % Re 

= ff/nn(6). (27) 

There is no net flux of material (liquid plus particles) through an arbitrary horizontal cross 
section of the suspension. Thus, if the particles settle, a net upward flux of fluid with velocity 
— (V)<&/(1 — $) results. The average relative velocity of particles to liquid is therefore 
u = (V)/(l — $). We now equate the expressions ( p6|) and (|27|), substitute u, and solve for 
(V): 

Here, we have replaced the liquid volume fraction e by the solid volume fraction $ = 1 — e. 
We obtain the dimensionless hindered settling function /hs(^) by division of this equation 
by V s , 

/hs($) = vs ~ 2 iur^y (29) 



For example, if we insert the functional form of / rm (e) found in Eq. (p3|) for the regime 
of Stokes flow and high dilution we obtain for the hindered settling function, 

/ hs ($) = (1 - $) 3 . (30) 

For the description of experimental data one often uses the phenomenological Richardson- 
Zaki correlation p2 |, 



/ hs ($) = (1 - $)", (31) 

where n takes values around 5. Theoretical arguments by Batchelor [[3(|[37]] lead to the 
expression 

/ hs ($) = 1-6.55$ (32) 

for small $. 

For comparison with the simulation results we have plotted relation ( |30"D in Fig. |5| as 
solid line. 
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IV. CONCLUSION AND OUTLOOK 



A. General advantages of a drag-force based approach 

The major advantage of a drag-force based algorithm as described here for the simulation 
of two-phase flows is the capability to simulate considerably large systems with relatively 
moderate computational requirements. At the same time one gets without efforts a physical 
behavior of the particulate phase circumventing the problems of continuum approaches with 
stress or the constitutive equation of the particulate phase. Such an algorithm is therefore 
a tool to obtain quick insight into collective effects of many particles. In fact, most of the 
known physical effects appear to be reproducible using only a very simple expression for the 
momentum exchange between phases. 



B. Specific problems of a drag- force based approach 

We list now some of the problems that we have encountered performing the present study 
that will be addressed in future research. 

Let us first consider the behavior of a single particle in an otherwise resting liquid, 
both initially at rest. If the particle density is higher then the density of the liquid, the 
particle will experience a net force in direction of gravity. As it falls, liquid will be dragged 
down with it in accordance with the momentum exchange rules and the requirements of the 
replacement of liquid as the local liquid fraction changes. As compared to the liquid velocity 
far away from the particle the local velocity is closer to the particle velocity. Consequently 
the isolated particle experiences a drag force which does not agree with the theoretical Stokes 
expression into which the particle velocity and the liquid velocity at infinity enter and the 
resulting terminal fall velocity of the particle is larger than the theoretical Stokes velocity.^ 

One way to address this problem is to alter the drag law (|T^|) by introducing a "drag- 
coefficient" different from 67r and dependent on the, albeit unphysical, ratio of particle size 
to the grid-spacing and physical parameters as the Reynolds number — this approach has 



been taken in ||38|| . However, the root of this problem lies in the treatment of the momentum 
exchange and indicates that a revision of its treatment is necessary. 

One could imagine first that a better treatment of the momentum exchange should 
proceed along the lines of approximating the liquid stress tensor on the particle surface, which 
may be possible in the low Reynolds number regime. However, the principal limitation of 
our ansatz is that the flow velocities in the vicinity of the particle surface are not accurately 
rendered as required for an estimation of the stress on the particle, since we account for the 
flow only in an average sense. 

Another possibility is to change the drag law. We have so far employed a simple Stokes 
expression. It is probably necessary to introduce a drag expression which additionally de- 
pends on the local particle density in addition to the local velocity difference between phases. 



2 This effect has been taken into account by using the measured terminal velocity of a single sphere 
in lieu of the theoretical Stokes velocity, for example in the scaling of the mean sedimentation 
velocities to obtain the hindered settling function. 
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We think that at least at high particle densities also particle rotation becomes important 
and must be included into the model. Moreover, we may have to use density dependent 
local viscosities in the "Navier- Stokes" equation. 

In conclusion, we should and will undertake an improvement of our algorithm to tackle 
some of the problems listed above in the future. However, even at the present state, we have 
a powerful tool at hand to assess cooperative effects of many-particle systems in which hy- 
drodynamic interactions play an important role. We have successfully modeled flow around 
impurities at small and moderate volume fractions and found reasonable agreement in the 
case of a sedimenting system. 
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TABLES 



L x system width 


3 


cm 


L y system height 


6 


cm 


z Stokes force reference length 


0.03 


cm 


pi liquid density 


1.09 


g/cm 3 


rj liquid shear viscosity 


0.4 


g/cm s 


f average particle radius 


0.015 


cm 


p polydispersity 


0.1 


1 


TABLE I. List of parameters for the simulation of flow through porous media. 


Deviating 


parameters are stated in the text. 






L x system width 


2.5 


cm 


L y system height 


5.0 


cm 


z Stokes force reference length 


0.05 


cm 


pi liquid density 


1.0 


/ 3 

g/cm 


rj liquid shear viscosity 


0.4 


g/cm s 


p p particle density 


2.53 


g/cm 3 


f average particle radius 


0.025 


cm 


p polydispersity 


0.01 


1 


k n repulsion constant 


2.5 x 10 4 


g/s 2 


7t dissipation constant 


80 


l/s 


p Coulomb friction factor 


0.3 


1 


g gravitational acceleration 


981 


cm/s 2 



TABLE II. List of parameters for the simulation of batch settling sedimentation. Deviating 
parameters are stated in the text. 
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FIGURES 



FIG. 1. 2D Fluid flow pattern through a random assembly of "spheres." In this figure the 2D 
area fraction is e = 0.15 and the particle based Reynolds number Re p = 0.01. Lines indicate the 
direction and the magnitude of the flow. Clearly visible are the effects of the mass conservation on 
the flow: The flow concentrates in regions with few particles and "engulfs" the particles on smaller 
scales. 

FIG. 2. We plot the product f p Re p of friction factor and particle Reynolds number vs. 
the particle Reynolds number Re p , at constant porosity e = 0.93 (= 1 — X)j(4/3)7rrf /L x L y z). 
Apart from small statistical fluctuations due to different initial placements of the particles, the 
curve is constant in the regime of Reynolds numbers <C 1, thus showing that the pressure drop is 
proportional to u as required by Darcy's law. 



FIG. 3. We plot the product f p Re p of friction factor and particle Reynolds number vs. the 
volume fraction e. Different curves correspond to different flow rates (Reynolds numbers) and 
to different values of the system "depth" z/f = 2 (+), z/f = 8/3 (□). A third set of runs has 
been made for the case of only frictional coupling between fluid and obstacles, setting the liquid 
volume fraction to one everywhere (O). The volume fraction e is calculated as 3D volume fraction, 
considering the particles as spheres and the simulation box as quadrilateral of extension L 

X X Ly X Z . 

FIG. 4. Typical particle configurations during a batch settling simulation at tVs/f = 5.2 
(a), 16 (b), 26 (c), 36 (d), 46 (e), and 56 (f) where t denotes time, Vs the Stokes velocity and f 
the average particle radius. The system size is chosen to be comparatively small L x /f = 80 and 
L y /L x = 2; 606 particles are visible and the 2D solid area fraction ^l i Trrf /L x L y is 0.10. Lines 
indicate the direction and amplitude of the fluid flow. Shades of gray denote the y-component of 
the particle velocity, dark particles move down and light ones up. It is interesting to see that there 
is substantial internal motion of particles in complicated vortex patterns. The particle settling is 
visible "on average" for example at the upper sedimentation front. 
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FIG. 5. Hindered settling function /h s = {V)/Vs for simulation series (i) (see text), including 
lubrication, void fraction and using periodic boundary conditions in x direction. The ordinate 
shows the effective 3D solid fraction <£, which is related to the 2D area fraction of particles in a 
simple way, given that the radius distribution of the particles is monodisperse: = (4/3)(f/z)$2D- 
The data points have been averaged over three runs with different initial conditions. The solid line 
is the form As = (1 — ^) 3 proposed in the text. 

FIG. 6. Ratios of sedimentation velocities vs. the effective 3D volume fraction <E> in the 
simulation computed under four different sets of simulation conditions (i)..(iv), for details see text. 
Diamonds (o) denotes the ratio the computed velocity disregarding lubrication (ii) to the velocity 
in the "full" simulation (i). Similarly, (+) denotes the ratio of series (hi) to (ii) — no lubrication 
and in (hi) additionally the liquid fraction set to 1. Boxes (□) denote the ratio (iv) to (i) for the 
difference between periodic and no-slip boundary conditions for the box walls. All data points have 
been computed using velocity averages over three runs with different initial conditions. 
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FIG. 1. 2D Fluid flow pattern through a random assembly of "spheres." In this figure the 2D area fraction is e = 0.15 and 
the particle based Reynolds number Re p = 0.01. Lines indicate the direction and the magnitude of the flow. Clearly visible are 
the effects of the mass conservation on the flow: The flow concentrates in regions with few particles and "engulfs" the particles 
on smaller scales. 
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FIG. 2. We plot the product f p Re p of friction factor and particle Reynolds number vs. the particle Reynolds number Re p , 
at constant porosity e = 0.93 (= 1 — ^\(4/3)-7rrf /L x L y z). Apart from small statistical fluctuations due to different initial 
placements of the particles, the curve is constant in the regime of Reynolds numbers <C 1, thus showing that the pressure drop 
is proportional to u as required by Darcy's law. 
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FIG. 3. We plot the product f p Re p of friction factor and particle Reynolds number vs. the volume fraction e. Different 
curves correspond to different flow rates (Reynolds numbers) and to different values of the system "depth" zjf = 2 (+), 
zjf = 8/3 (□). A third set of runs has been made for the case of only frictional coupling between fluid and obstacles, setting 
the liquid volume fraction to one everywhere (<>). The volume fraction e is calculated as 3D volume fraction, considering the 
particles as spheres and the simulation box as quadrilateral of extension L X X Ly X Z. 



FIG. 4. Typical particle configurations during a batch settling simulation at tVs/f = 5.2 (a), 16 (b), 26 (c), 36 (d), 46 (e), 
and 56 (f) where t denotes time, Vs the Stokes velocity and f the average particle radius. The system size is chosen to be 
comparatively small L x /f = 80 and L y /L x = 2; 606 particles are visible and the 2D solid area fraction 7rrf /L x L y is 0.10. 
Lines indicate the direction and amplitude of the fluid flow. Shades of gray denote the y-component of the particle velocity, 
dark particles move down and light ones up. It is interesting to see that there is substantial internal motion of particles in 
complicated vortex patterns. The particle settling is visible "on average" for example at the upper sedimentation front. 




FIG. 5. Hindered settling function /h s = (V) /Vs for simulation series (i) (see text), including lubrication, void fraction and 
using periodic boundary conditions in x direction. The ordinate shows the effective 3D solid fraction which is related to the 2D 
area fraction of particles in a simple way, given that the radius distribution of the particles is monodisperse: $ = (4/3)(f/z)$2D- 
The data points have been averaged over three runs with different initial conditions. The solid line is the form /h s = (1 — $) 3 
proposed in the text. 



1.2 
1.15 

1.1 
1.05 
1 



0.95 



0.9 
0.85 
0.8 

0.05 0.1 0.15 0.2 

3> 




FIG. 6. Ratios of sedimentation velocities vs. the effective 3D volume fraction $ in the simulation computed under four 
different sets of simulation conditions (i)..(iv), for details see text. Diamonds (o) denotes the ratio the computed velocity 
disregarding lubrication (ii) to the velocity in the "full" simulation (i). Similarly, (+) denotes the ratio of series (iii) to (ii) 
— no lubrication and in (iii) additionally the liquid fraction set to 1. Boxes (□) denote the ratio (iv) to (i) for the difference 
between periodic and no-slip boundary conditions for the box walls. All data points have been computed using velocity averages 
over three runs with different initial conditions. 



